Structural, elastic, electronic, optical and vibrational properties of single-layered, bilayered and bulk molybdenite MoS2-2H

An extensive density functional theory analysis is reported of the crystal-chemical properties, electronic band structure, and optical (absorption and electron energy-loss spectra) and phonon properties of bulk and mono- and bilayer molybdenite, a transition metal dichalcogenide with important applications in semiconductor technology. This mineral is characterized by virtually infinite bidimensional layers of MoS2 (with covalent Mo—S bonds) held together by weak dispersive forces, explaining the easy cleavage of the (001) surfaces.


Introduction
Molybdenite (MoS 2 ) is a sulfide mineral that crystallizes in the hexagonal crystal system (space group P6 3 /mmc) under standard pressure and temperature conditions. The unit cell of molybdenite is shown graphically in Fig. 1, presenting two formula units (Z = 2) of molybdenum disulfide, each one forming a layer where the atoms are linked by mixed covalent/ ionic bonds. The molybdenum atom is sixfold coordinated, with a trigonal prism arrangement of the S 2À ions. These layers are extended indefinitely along the a and b axes and stacked along the c-axis direction. The whole structure is held along the [001] direction by weak long-range (van der Waals) interactions. Because the crystal structure is made up of stacked layers, molybdenite is subject to polytypism, and the structural model in Fig. 1 represents the MoS 2 -2H polymorph stable under standard conditions of pressure and temperature (1 atm, 298.15 K; 1 atm = 101 325 Pa).
Molybdenite is a widely employed mineral phase in several and manifold geological, industrial and technological applications. It is the main ore mineral for the extraction of metallic molybdenum (Hess, 1924) and an element used in alloy steels, superalloys and pigments, just to cite a few examples (Kropschot, 2010). Bulk molybdenite can be used as a dry lubricant because of its very low friction coefficient (Martin et al., 1993(Martin et al., , 1994, even under high pressure (Wu et al., 2022).
In addition, molybdenite was the first discovered semiconductor material in the early 20th century (Aminoff & Broome, 1935), and it has drawn much attention in the past decade for its use as a two-dimensional (2D) material in fieldemission transistors (FETs) and other electronic devices, such as solar cells and light-emitting diodes (Radisavljevic et al., 2011;Fontana et al., 2013;Sebastian et al., 2021;Wang et al., 2021). Thanks to the weak interactions along the [001] direction, molybdenite is easily cleavable, and a monolayer of molybdenite was exfoliated by Joensen et al. (1986) using intercalated lithium. It was shown that a single layer of MoS 2 has superior properties to graphene for the fabrication of FETs and optoelectronic devices because it presents a gap in the band structure (Radisavljevic et al., 2011). This band gap is tunable by increasing or decreasing the number of layers of the material (Yang & Li, 2020), and even changes from an indirect gap of about 1.23 eV in bulk MoS 2 -2H (Kam & Parkinson, 1982) to a direct gap of about 1.8 eV in the monolayer (Kuc et al., 2011). This indirect-direct band-gap transition is fundamental to enhancing light emission and optical absorption, and to inducing photoluminescence in the monolayer of molybdenite (Mak & Shan, 2016).
Several experimental and theoretical reports have investigated the structural, electronic and vibrational properties of bulk and layered molybdenite to provide a better understanding of the quantum confinement effects originating from 'layering' the mineral and to drive the development of MoS 2based materials for photonics and optoelectronic applications. However, most of them focused only on specific features, such as the band gap (Kuc et al., 2011;Kam & Parkinson, 1982), the dielectric properties (Beal & Hughes, 1979;Kumar & Ahluwalia, 2012), the elasticity (Feldman, 1976;Peelaers & Van de Walle, 2014) and the vibrational properties Funke et al., 2016;Wieting & Verble, 1971). Some of these reports also showed conflicting results, such as a shift of some Raman modes when MoS 2 is prepared in layers.
Because of the relevance of this material in the next generation of devices, and to provide a detailed and comprehensive cross-correlated analysis of the different properties of molybdenite and its layered derivatives, we have performed a series of ab initio simulations at the density functional theory (DFT) level of the structures of the MoS 2 -2H bulk, the monolayer (here labelled as MoS 2 -1L) and the bilayer (MoS 2 -2L). This choice is dictated by previous theoretical analyses (Kuc et al., 2011), which highlighted the occurrence of the direct band gap only in the monolayer, with bilayer MoS 2 behaving at the electronic level similarly to the bulk mineral. In the following, an analysis of the structures and elastic, electronic, dielectric and phonon properties is provided, cross-correlated, and discussed against previous results reported in the scientific literature. Most important in our work is that all the simulations were performed with a correction to include the van der Waals interactions in the physical treatment, a non-trivial choice when dealing with layered materials such as molybdenite, other transition metal dichalcogenides and phyllosilicates (Moro et al., 2016).

Computational methods
The present simulations have been conducted by means of the Vienna ab initio Software Package (VASP) code (Kresse & Furthmü ller, 1996;Kresse & Hafner, 1993) using the density functional PBE (Perdew et al., 1996). The projector augmented-wave (PAW) basis set was expanded using a cutoff of 600 eV throughout the different calculations (Kresse & Joubert, 1999). The contribution of the weak van der Waals forces was included according to the DFT-D3 approach (Grimme et al., 2011), where the first two summations are over all N atoms in the cell and the third one is over the translations of the unit cell L = (l 1 , l 2 , l 3 ). C 6ij and C 8ij are the dispersion coefficients of the atom pair ij, and r ij,L is the distance between atom i in the reference cell (L = 0) and atom j in the cell L. The function f d,n is the Becke-Jonson damping function, which is given by the expression with s 6 = 1 and a 1 , a 2 and s 8 as variable parameters depending on the geometry of the system. Geometry optimization of the bulk (MoS 2 -2H) was conducted by means of an iterative procedure involving the successive minimization of the a lattice parameter and c/a ratio, starting from the experimentally refined unit cell (Schö nfeld et al., 1983). Within this procedure, the tolerance on forces was set at 10 À4 eV Å À1 . The molybdenite monolayer (MoS 2 -1L) and bilayer (MoS 2 -2L) structures were modelled as slabs containing just one and two structural layers, respectively, placed in a unit cell with the a axis equal to the refined a parameter of the bulk MoS 2 -2H. A vacuum region along the [001] direction of 20 and 30 Å thickness for MoS 2 -1L and MoS 2 -2L, respectively, was employed to avoid interaction between replicas of layers in neighbouring cells along the c axis. The sampling in the first Brillouin zone for the bulk mineral was conducted on a Monkhorst-Pack À-centred grid (Monkhorst & Pack, 1976) with 17 Â 17 Â 9 k-points for the geometry optimization procedure, whereas a finer mesh of 33 Â 33 Â 9 k-points was employed for the calculation of the electronic band structure, density of states and dielectric function. For the layered systems MoS 2 -1L and MoS 2 -2L, due to the increased c lattice parameter, the grids were reduced along the [001] direction to just one k-point, leading to 33 Â 33 Â 1 and 33 Â 33 Â 1 k-point meshes. For all the calculations, the threshold controlling the convergence on the total energy between two consecutive self-consistent field steps was set to 10 À8 eV. Electronic properties (band structure and density of states) at the GW level of theory were obtained through generating maximally localized Wannier functions, which were calculated using the Wannier90 code (Pizzi et al., 2020). The frequency-dependent dielectric function was calculated with zero momentum transfer [q = 0; q = (4/ )sin, where is half the scattering angle and is the wavelength of the incident radiation] using both the independent particle approximation (without local field effects) (Gajdoš et al., 2006) and the random phase approximation, the latter performed on top of a single-shot G 0 W 0 calculation (Shishkin & Kresse, , 2007. The frequency-dependent dielectric function was then compared with known measures and models reported in the scientific literature. Phonon properties were calculated using the Phonopy package (Togo & Tanaka, 2015), calculating the À-point (q = 0) frequency and the phonon dispersion relations (q 6 ¼ 0) using finite displacements and the modified Parlinski-Li-Kawazoe method (Parlinski et al., 1997) on 3 Â 3 Â 3 and 3 Â 3 Â 1 supercells for the bulk and layered molybdenite models, respectively.

Results and discussions
3.1. Molybdenite structure Geometry optimization is the first step necessary prior to any further analysis to assess the quality of the simulation approach. The structural results for the different molybdenite models (bulk, monolayer and bilayer) are reported in Table 1, alongside previous experimental and theoretical determinations. For the bulk mineral, our simulation at the PBE-D3 level of theory is correctly able to reproduce the structural features observed from the X-ray diffraction (XRD) refinements of Bronsema et al. (1986). In detail, we observed a small underestimation of the a and c lattice parameters of about À0.35 and À1.74%, respectively, which is an expected result at 0 K. The bond lengths are in good agreement with the XRD results, with a slight underestimation of about À1.88 and À0.46% for the S-S and Mo-S bonds, respectively. The structural data are also in good agreement with the experimental lattice parameters obtained via scanning tunnelling microscopy, a = 3.160 Å and c = 12.294 Å (Schö nfeld et al., 1983). Furthermore, our theoretical Mo-S and S-S bond lengths (2.399 and 3.149 Å , respectively, at the DFT-D3 level) are in good agreement with those measured experimentally (2.41 and 3.19 Å , respectively;Schö nfeld et al., 1983).
The binding energy between the MoS 2 layers was calculated according to the following formula: where E bulk is the total energy of the bulk mineral and E 1L is the energy of a single layer (the factor 2 considers the presence of two layers in the bulk unit cell). Our theoretical value (ÁE bind = 47 meV per atom = 0.456 J m À2 ) is in very good agreement with both the ÁE bind = 0.52 (4) J m À2 obtained by Weiss & Phillips (1976) from the experimental specific surface energy of molybdenite and the binding energy value of 0.55 (13) J m À2 measured by Fang et al. (2020) using an in situ peeling-to-fracture method. In the latter study, the authors performed a theoretical simulation of the mechanical exfoliation of molybdenite with the VASP code using the optB86b-vdw kernel, which is a modified version of the vdW-DF approach. The results provided ÁE bind = 0.422 J m À2 , which is slightly lower than the value obtained from our simulations. We infer that this small deviation is due to the different computational settings employed by Fang et al. (2020), especially the different description of the van der Waals long-range interactions. We also calculated the cohesive energy of bulk molybdenite using the formula and of the layered structures using research papers where E Mo A and E S A are the atomic energy of the isolated Mo and S atoms, respectively, n is the number of layers, and E nL is the energy of the layered system. The cohesive energy (5.438 eV per atom) is in good agreement with the experimental value of 5.18 eV per atom reported by Raybaud et al. (1997). These results are also in line with those of Bučko et al. (2013), who obtained ÁE coh = 5.37 eV per atom. As explained by those authors, this slight overestimation is mainly related to a high contribution from the correction for the dispersive forces, since the sole PBE functional provided a cohesive energy closer to the experimental value. The ÁE coh value that we obtained from the DFT-D3 approach seems to confirm this hypothesis.
Comparison with the results using the standard PBE functional  indicates that the inclusion of dispersive forces, at least via a posteriori corrections, plays a fundamental role in determining the physical-chemical properties of MoS 2 -2H and other layered materials, in agreement with previous statements on the topic (Cutini et al., 2020;Ulian et al., 2013Ulian et al., , 2021. Indeed, the lack of an appropriate treatment results in an overly overestimated c lattice parameter and extremely low bulk modulus, as shown by . The uncorrected PW91 functional provided a good description of the structural information directly related to the MoS 2 layers, i.e. the a lattice parameter, but the interlayer binding energy dropped to the negligible value of 3 meV per atom (À94%), resulting in a dramatic increase of the c lattice parameter (25% with respect to the PW91-D2 results). Hence, since most of a mineral's properties (optical, electronic and so on) are strictly related to its structure, it is of the utmost importance to use a theoretical approach that considers all the relevant physical forces.
One of the first attempts to include dispersive forces in the physical description of layered structures (e.g. graphite and molybdenite) was provided by Rydberg et al. (2003), who included these kinds of forces in a nonlocal correlation density functional called vdW-DF, coded in the PWscf program. The structural results (lattice parameters and bulk modulus) provided were in reasonably good agreement with the available experimental data Bronsema et al., 1986). Compared with the results of the present study, the vdW-DF approach provided a slightly overestimated unit-cell volume due to larger a and c parameters. More recently, in the work of Bučko et al. (2013), who employed the VASP code and a DFT-D2 correction for the dispersive forces, the calculated lattice parameters of the unit cell were a = 3.19 Å and c = 12.42 Å , which are slightly larger than those obtained in the present study, in particular for the c axis. This discrepancy is mainly related to the different correction for the dispersive forces, because the original D2 scheme of Grimme (2006) is more empirical than the D3, where the C 6 parameters are adjusted for each atom according to the local chemical environment (charge density) in which it is located. Also, the 8 Â 8 Â 8 k-point mesh employed by Bučko et al. (2013) has fewer sampling points than the one adopted here (17 Â 17 Â 9), although the authors considered a very large kinetic energy cut-off (i.e. quality of the basis set) of 1500 eV.
Our simulation results for bulk MoS 2 -2H are also in good agreement with previous theoretical results at the DFT level reported by . In that work, the authors employed the PWscf package in QuantumEspresso (Giannozzi et al., 2009), plane-wave basis sets with ultrasoft pseudopotentials (Vanderbilt, 1990) and the PW91 density functional (Perdew et al., 1992), corrected with the DFT-D2 scheme. In particular, the unit-cell lattice and internal geometry of bulk molybdenite obtained with our approach were within about  Table 1 Calculated lattice parameters a and c (Å ), S-S and Mo-S bond distances (Å ), S-Mo-S bond angle Â ( ), interlayer distance d I (Å ), internal parameter z (dimensionless), bulk modulus (K 0 , GPa), interlayer binding energy ÁE bind (meV per atom) and cohesive energy ÁE coh (eV per atom) for bulk (MoS 2 -2H) and layered (MoS 2 -1L and MoS 2 -2L) structures, compared with previous theoretical and experimental (Exp.) results.
The asterisk (*) marks fixed parameters in the present theoretical simulations.

System
Method 1% of those previously simulated, and no difference in the binding energy ÁE bind was observed. We note that the formation energy from our simulation is about 6% more than that calculated at the PW91-D2 level, which may be due to the different choice of the density functional and basis sets. Compared with the bulk mineral, the molybdenite monolayer (MoS 2 -1L) and bilayer (MoS 2 -2L) structures showed very small variations regarding their internal geometry, with differences not higher than 0.1%. This is a typical trend in layered minerals and in materials whose structural features do not vary significantly when simulated in layers (Ulian et al., 2018). This can also be evinced from the cohesive energy of both the mono-and bilayer models, which were within 1% of that of bulk MoS 2 -2H. It is interesting that the difference between the cohesive energy of MoS 2 -1L and MoS 2 -2H corresponds to the binding energy per MoS 2 unit, which is twice the ÁE bind per atom. The same applies to the bilayer slab of molybdenite. Finally, and as expected, the binding energy between the two layers of the molybdenite bilayer structure is 22 meV per atom = 0.213 J m À2 , which is half the ÁE bind value of the bulk mineral because in the latter each MoS 2 sheet interacts with one above and one below, whereas there is only a single layer-to-layer interaction in the bilayer model. The above observations and discussion are in line with the results reported by , who compared the structural properties of bulk and monolayer molybdenite. However, those authors did not consider a possible bilayer structure in their simulations.

Elasticity
To provide a further assessment of the quality of our simulation approach, we calculated the second-order elastic moduli of bulk molybdenite. The elasticity of bulk molybdenite was obtained using a finite strain approach, calculating the stiffness moduli from the stress-strain relationship (Yu et al., 2010;Nye, 1957). For a hexagonal structure there are five independent elastic moduli that, according to the 6 Â 6 matrix notation of Voigt, can be expressed as where C 66 = (C 11 À C 12 ) /2 and the dots are null moduli. Three strain patterns were necessary to obtain the elastic constants and, for each of them, five equally spaced strain magnitudes between À0.015 and 0.015 were employed. For the calculation of the elastic moduli, the unit cell vectors a and c were oriented parallel to the x and z Cartesian axes, respectively, which is the standard crystallographic orientation of a hexagonal unit cell proposed by the Institute of Radio Engineering (Brainerd, 1949). The simulation results are presented in Table 2 and compared with previous experimental and theoretical data. Note the strong anisotropy between the C 11 = 231.82 GPa and C 33 = 63.47 GPa stiffness matrix components, which are related to elastic deformations along the a and c axes, respectively. As previously mentioned, this is a consequence of the different bonding scheme in molybdenite, where the Mo and S atoms are linked by strong covalent bonds and the MoS 2 layers, stacked along the [001] direction, are held together by van der Waals interactions. A similar anisotropic behaviour has been observed for other layered phases as well (Gatta et al., 2015;Ulian et al., 2014). These theoretical findings are in line with the approximate values calculated from different experimental measurements (neutron dispersion curves and compressibility data from XRD refinements) by Feldman (1976). Our results show an underestimation of about À2.6% for the C 11 modulus and an overestimation of the C 13 stiffness component of about 22%, but they nevertheless fall within the large uncertainties associated with the experimental data reported and discussed by Feldman (1976). In addition, we are aware that the C 12 stiffness component is positive in our work and negative in the cited study, a discrepancy caused by the different orientation of the Cartesian axes relative to the unit cell of the mineral (see above).

Electronic properties
The electronic band structure and density of states for molybdenite bulk (MoS 2 -2H), calculated along the K-G-M-K-H-L-A-H path in the first Brillouin zone, are reported in Fig. 2, and the same results are presented in Fig. 3 for the MoS 2 -1L (monolayer) and MoS 2 -2L (bilayer) structures. It is possible to note from these results that there are four groups of bands in the calculated structure and density of states. The first group is related to valence bands located at low energy, between À12 and À15 eV, related mainly to the 3s orbital of the sulfur atom. The second group of valence bands, below the Fermi energy and separated from the first one by a large gap of about 7 eV, is the result of a visible hybridization of the 3p and 4d orbitals of S and Mo, respectively, with a small contribution from the 5s orbital of molybdenum around À5 eV. The third group is above the Fermi energy, representing the first range of conduction bands of molybdenite that are predominantly the result of the 4d orbitals of Mo and, to lesser extent, the 3p orbitals of S. The last group of (conduction) bands is located above ca 5 eV, given by a mixed contribution (hybridization) of the 5s and 5d states of  molybdenum and the 3p orbitals of sulfur. Since the d orbitals are the main contributor to the bands located near the band gap, the bands are almost flat. While these features are shared between the different structures (bulk and monolayer), there are important differences arising from the thickness of the materials. For instance, in the bulk MoS 2 -2H, the maximum energy of the highest occupied valence band is located on the À point and the minimum energy of the lowest occupied conduction band is on a point in the K ! À path. This means that bulk molybdenite is a semiconductor with an indirect band gap of 0.79 eV. The band gap becomes a wider direct one, E g = 1.90 eV, located on     the K point for a monolayer of the mineral. However, even the bilayer MoS 2 structure reverts the band gap from the direct K-K 0 to the indirect one observed for the bulk mineral, albeit the gap is slightly higher (1.22 eV).
The present results are in very good agreement with previous theoretical ones reported in the literature. For example, Kumar & Ahluwalia (2012), using Troullier-Martin norm-conserving pseudopotentials as coded in the SIESTA software, calculated a band gap for bulk molybdenite of 0.75 eV by means of the local density approximation (LDA) functional, which increased to 1.05 eV when the authors employed the DFT functional of Perdew et al. (1996). For the MoS 2 -1L monolayer the obtained values were higher, namely 1.89 and 1.55 eV at the LDA and PBE levels of theory, respectively.  performed similar simulations by means of the projector augmented-wave basis set, calculating the band gap of bulk MoS 2 using LDA (or the generalized gradient approximation) as 0.72 eV (0.85 eV), whereas they obtained 1.87 eV (1.58 eV) for the molybdenite monolayer.
However, as also reported by Kumar & Ahluwalia (2012), all theoretical simulations on bulk MoS 2 -2H underestimated the band gap, as the experimental investigation by photocurrent spectroscopy assessed this value at 1.23 eV (Kam & Parkinson, 1982). This is a common issue related to both LDA and GGA (generalized gradient approximation) DFT functionals. To overcome this drawback, in the present work the GW approximation was employed, which allows us to describe the quasiparticle electronic properties. 4(e)] and the MoS 2 -2L bilayer [Figs. 4(c) and 4( f )]. The calculated band gaps are 1.14 eV (indirect), 2.57 eV (direct) and 1.88 eV (indirect) for the bulk mineral, monolayer and bilayer MoS 2 , respectively. Hence, the GW approach provides a better description of the band gaps of the MoS 2 -2H bulk and the layered structures. The results for the monolayer are also in very good agreement with the theoretical results of Zibouche et al. (2021), who calculated the band gap (2.68 eV) using the Sternheimer equation .

Optical response function
The optical properties are due to the electronic transition from the occupied to the unoccupied state, in other words they are two-particle excitations. At the DFT level, the key quantity that is calculated is the frequency-dependent complex dielectric function "(!) = " 1 (!) + i" 2 (!), with " 1 and " 2 the real and imaginary parts, respectively. This knowledge is very useful because the real component of "(!) provides information on the transmission of electromagnetic waves through the medium, whereas the " 2 (!) component is related to interband electronic transitions, i.e. single-particle excitations (Raether, 1980). The imaginary part of " is calculated as a summation over empty states, according to the equations reported by Gajdoš et al. (2006), whereas " 1 (!) was obtained using the Kramers-Kronig relation. In addition, the dielectric function can also be calculated with the electric vector oscillating either parallel or perpendicular to the c axis, because each considered phase belongs to the hexagonal system. This is calculated with the following formulae: and Other properties of interest can be derived from the complex dielectric function, such as the refractive index n(!), the extinction coefficient (!), the absorption coefficient (!), the energy-loss function EELS(!) and the reflectivity R(!), according to These optical properties were evaluated for each molybdenite model for the MoS 2 -2H, MoS 2 -1L and MoS 2 -2L structures, considering a high number of unoccupied bands to increase the accuracy of the results. Fig. 5 reports the dielectric function calculated in the photon energy range 0-30 eV, subdivided into real (" 1 ) and imaginary (" 2 ) parts. Bulk molybdenite [Figs. 5(a) and 5(b)] shows excellent agreement with the data measured by Beal & Hughes (1979)   the electric field oscillating perpendicular to the (001) plane. In this case, the calculated optical properties were not shifted to give a better match to the experimental findings, as was done in previous work (Kumar & Ahluwalia, 2012). Three main peaks can be observed in the low-energy region (1-5 eV) of the imaginary part " 2 , labelled A (2.76 eV), B (4.21 eV) and C (5.31 eV) according to the nomenclature proposed by Kumar & Ahluwalia (2012). The A peak corresponds to the C exciton (Funke et al., 2016). There is good agreement with the theoretical data reported by Kumar and Ahluwalia, who calculated A = 2.9 eV, B = 4.5 eV and C = 5.1 eV, but our computational setting allows us to detect the exciton at 1.72 eV, corresponding to the K-K 0 transition from the topmost valence (v1) band to the first conduction band (c1), marked with a black arrow in Fig. 5(a). However, the intensity of the signal is lower than the experimental one, and the second expected excitonic transition c2-v1 is not visible because of the overlap with the A peak. The present data are also in line with thermal ellipsometry data at 35 K measured by Le et al. (2019), who measured A = 2.91 eV and the K-K 0 transition at 1.96 eV.
MoS 2 -1L presents four peaks in the imaginary part of the complex dielectric function, labelled A (2.82 eV), B (3.70 eV), C (4.39 eV) and D (5.54 eV) in Fig. 5(c), and they are in agreement with the results of Kumar & Ahluwalia (2012), i.e. A = 2.9 eV, B = 3.8 eV, C = 4.5 eV and D = 5.5 eV. In this case, the B peak is a new feature that appears because of the singlelayer structure, whereas the other signals (A, C and D) are slightly shifted at higher energy with respect to the A, B and C peaks of the mineral bulk. A single excitonic K-K 0 transition is also clearly visible at about 1.94 eV, in excellent agreement with that found at 1.9 eV by Funke et al. (2016), which was derived from spectroscopic ellipsometry measurements. A second peak at 2.05 eV was also measured by the same authors, which is related to spin-orbit coupling, i.e. a splitting of the topmost valence band in the single-layered MoS 2 unit at about 140-150 meV. An example of this effect on the bands of 2D molybdenite can be seen in the experimental and theoretical work of Song et al. (2019) and Moynihan et al. (2020). Since spin-orbit coupling was not included in the present simulations, this second peak (K-K 0 transition) is not visible in the dielectric function. Funke et al. (2016) found at about 3 eV the excitonic transition here labelled as A, in line with our calculations.
When the molybdenite structure is made of two layers, the shape of the dielectric function reverts to a form similar to that of bulk MoS 2 -2H [see Fig. 5(e)], presenting only three major peaks, A (2.73 eV), B (4.32 eV) and C (5.36 eV), in the range 0-5 eV.
In general, all the observed optical transitions occur in an energy range corresponding to electronic transitions between the p valence (occupied) bands of sulfur and the d conduction (unoccupied) bands of molybdenum (see Fig. 3), as also observed in previous work (Beal & Hughes, 1979;Kumar & Ahluwalia, 2012). According to the literature, the A peak in Figs. 5(a), 5(c) and 5(d) is mainly due to nearly parallel bands in the M ! À path, i.e. band nesting. The low-energy excitons are instead related to electronic transitions from the Mo d z orbital, which is split at the K point in the Brillouin zone because of spin-orbit interaction.
The real part of the dielectric function of the MoS 2 -2H structure is also in very good agreement with the experimental findings (Beal & Hughes, 1979;Funke et al., 2016), with values at zero energy of " k 1 = 10.42 eV and " ? 1 = 16.43 eV that are in line with those measured from experiments, e.g. " k 1 = 10.42 eV and " ? 1 = 16.8 eV (Beal & Hughes, 1979). These values are reduced when going from bulk to the monolayer (" k 1 = 3.29 eV and " ? 1 = 5.49 eV) and the bilayer (" k 1 = 4.53 eV and " ? 1 = 7.25 eV), showing a trend that increases with the number of MoS 2 units in the structure. Our results are also in line with the theoretical simulations of Kumar & Ahluwalia (2012), who calculated " k 1 = 8.9 eV and " ? 1 = 12.8 eV for the MoS2-2H bulk, and " k 1 = 3.0 eV and " ? 1 = 4.8 eV for the molybdenite monolayer. A better agreement can be observed by comparing our results with the theoretical ones of Ben Amara et al. (2016), who obtained " k 1 = 8.3 eV and " ? 1 = 15.4 eV from PBE simulations.
The calculated electron energy-loss spectroscopy (EELS) spectra (Fig. 6) provide further information on the plasmon modes, i.e. the collective oscillation of valence or conduction electrons in a material. These modes are related to transitions from the and bonding states to the respective antibonding states * and *, which occur in the EELS spectra as lowenergy and high-energy + plasmons (interband plasmons). Such plasmon modes are present in the EELS spectrum when the real part of the dielectric function " 1 (!) crosses zero with a positive slope. In bulk molybdenite (MoS 2 -2H), the plasmon falls at 8.52 eV, whereas the + one can be seen at 22.55 eV for E ? c. For an electric field oscillating parallel to the c axis, the + plasmon presents two peaks at 22.98 and 23.17 eV, and no other well defined signals are visible in the spectrum. A similar trend is observed for layered MoS 2 , but both the monolayer and bilayer systems show a red shift of these plasmon signals (for E ? c), at 8.06 eV () and 16.17 eV (+) for MoS 2 -1L, and 8.31 eV () and 17.44 eV (+) for the MoS 2 -2L model. For the electric field parallel to the c axis, the + plasmon falls at 15.55 and 17.01 eV for MoS 2 -1L and MoS 2 -2L, respectively. These theoretical results are in very good agreement with the experimental EELS measurements conducted via scanning transmission electron microscopy (STEM) by Moynihan et al. (2020), who obtained for MoS 2 -2H 8.6 and 23 eV for the and + plasmons, respectively. Compared with the previous DFT simulations at the PBE level (Kumar & Ahluwalia, 2012), our results are in general agreement, with the plasmon falling at lower energies than previously reported, i.e. 9.2 and 8.6 eV for the bulk and monolayer structures, respectively.
The optical absorption is also reported in Fig. 6, in the photon energy range 0-30 eV. The simulation results for bulk molybdenite [ Fig. 6(b)] are in excellent agreement with the experimental ones of Beal & Hughes (1979) up to about 15 eV, and then there is some deviation between the curves due to the data extrapolation performed by those authors. In the region of the visible spectrum [Vis, 380-780 nm, see the inset

Phonon properties
Bulk molybdenite (space group P6 3 /mmc, point group D 4 6h ) has six atoms in the unit cell, resulting in 18 degrees of freedom that, in the À point and from group-theory analysis, have the following irreducible representations (irreps): 2E 2g E 1u E 1g being the acoustic and optical modes, respectively. Modes with character A or B are non-degenerate, whereas the E modes are doubly degenerate. Vibrational motions that are symmetric (gerade, g) and anti-symmetric (ungerade, u) with respect to the inversion centre of the crystal are active in Raman and infrared spectroscopies, respectively.
The single layer of MoS 2 (P6m2 space group, D 3h point group) has only nine degrees of freedom and different zonecentre irreps, with À a = A 00 2 E 0 and À o = A 0 1 A 00 2 E 0 E 00 . If we consider a correlation with the D 6h point group, it is possible to associate the E 00 and A 0 1 modes of the monolayer with the E 1g and A 1g modes of the bulk, respectively. Instead, while the bilayer MoS 2 -2L has the same degrees of freedom as the mineral bulk, the P3m1 space group (D 3 3d point group) leads to different irreps,  . Red and blue lines refer to the on-plane (xx, E ? c) and out-of-plane (xx, E || c) electric field oscillations, respectively. Experimental data from Beal & Hughes (1979) are shown for direct comparison. The green dashed boxes in panels (b), (d) and ( f ) show the part of the absorption spectrum in the visible region, which is reported in the upper right inset.
where À a = A 2u E u and À o = 3A 1g 2A 2u 2E u 3E g . Again, gerade (A 1g and E g ) and ungerade modes (A 2u and E u ) are active in Raman and IR, respectively.
À-point frequencies are reported in Table 3, alongside selected results from both experiments (Funke et al., 2016;Wieting & Verble, 1971) and theoretical simulations Jimé nez Sandoval et al., 1991). The phonon dispersion curves for each model along the K-À-M-K-H path in the first Brillouin zone are shown in Fig. 7. No negative value of the phonon branch was found, meaning that the bulk and the free-standing layered structures are stable. In general, there is very good agreement between the different results, with mean absolute differences less than 3 cm À1 . Two Raman modes are quite sensitive to the number of layers in the structure, i.e. the No. 7 E 2g and No. 10 A 1g modes in MoS 2 -2H, the former decreasing and the latter increasing with the number of layers of the 2D material, in agreement with previous experimental observations (Funke et    thickness of the molybdenite. This is due to the stacking of the MoS 2 layers along the c axis, which affects the intralayer bonding (see Table 1) and hence the force constants of the vibrational motion of the atoms. Interestingly, according to Wieting & Verble (1971), the frequency difference between the E 1u and E 2g modes in the bulk is associated with the interlayer interaction, which is quite low in our simulations (no more than 3 cm À1 ). This extends to the bilayer model MoS 2 -2L, considering the E u and E g vibrations. However, while the absolute frequency shift of the A 1g ! A 0 1 mode is quite in line with the experimental one (ca 4 cm À1 ) of Lee et al. (2010), the E 2g ! E 0 mode is much less affected by the thickness of the material, in agreement with previous theoretical observations . We suppose that this is due to the absence of any interaction of the MoS 2 monolayer or bilayer with the support (e.g. SiO 2 or sapphire) that is commonly employed in experimental measurements.
Regarding the phonon dispersion (Fig. 7), there is a satisfactory agreement with previous inelastic neutron scattering data measured by Wakabayashi et al. (1975) along the À-M path ([100] direction) for MoS 2 -2H. No significant differences in the phonon density of states were observed in our simulations by changing the number of layers in the system.
Quite good agreement was found by comparing the present results with those of , who performed PW91-D simulations for the bulk and monolayer of molybdenite for both the phonon branches and the density of states. Their GGA+D approach resulted in E 1g = 286.6 cm À1 , E 2g = 378.5 cm À1 and A 1g = 400.2 cm À1 for MoS 2 -2H, and E 0 = 380.2 cm À1 and A 0 1 = 406.1 cm À1 for MoS 2 -1L. It is worth noting the opposite behaviour of the A 1g ! A 0 1 frequency shift, which is explained by the different unit-cell parameters between theory and experiment.

Conclusions
In this work, with the density functional theory method, using plane-wave basis sets, the PBE functional and the DFT-D3 scheme to include van der Waals interactions, we have investigated different properties of bulk molybdenite MoS 2 -2H and its monolayer (MoS 2 -1L) and bilayer (MoS 2 -2L) structures. The scope was to provide a detailed cross-correlated analysis of these materials on the atomic scale, both to increase the knowledge of them and as a comparison basis for future work.
We have noted subtle variations in the crystal chemistry of MoS 2 when simulated in layers, in particular the bond distances and angles, with binding energy values in agreement with experimental measurements that employed the peelingto-fracture method. The selected PBE-D3 approach also described the structure of bulk MoS 2 -2H, with the expected unit-cell volume at 0 K smaller than that from roomtemperature XRD refinements. The stiffness of bulk MoS 2 -2H was also in line with the only experimental work reported in the literature, with a slight overestimation of the C 33 elastic modulus and an underestimation of the C 13 one as a result of the approximations introduced in the simulations, chiefly the lack of thermal effects (athermal results at 0 K without any zero-point contribution).
The band gaps of bulk molybdenite, MoS 2 -1L and MoS 2 -2L calculated at the GGA (GW) level were an indirect 0.79 eV (1.14 eV), a direct 1.90 eV (2.57 eV) and an indirect 1.22 eV (1.88 eV), respectively, in line with both experimental data and theoretical simulations that considered quasi-particle excitations. We observed from the atom-projected density of states that the 3p orbitals of S hybridize with the 4d ones of Mo at the top of the valence band and at the bottom of the conduction band, whereas the core states are due to the 3s orbitals of the sulfide ions. We expect that hybrid functionals, such as HSE06, could perform similarly but with much higher computational costs, especially when using plane-wave basis sets.
The calculated complex dielectric function showed the direct (zero-momentum q) electronic transitions (imaginary part) and the plasmonic resonances (EELS spectrum) in the photon energy range 0-30 eV, which agree with recent experimental measurements. Different absorption peaks in the visible portion of the light spectrum were observed according to the number of layers of molybdenite.
Finally, we studied the IR and Raman vibrations at the zone centre and the phonon dispersion relations of the three MoS 2 models (bulk, monolayer and bilayer). Our simulations confirm the experimental evidence that the No. 7 E 2g mode decreases and the No. 10 A 1g mode increases with the number of layers of the bidimensional material. However, the extent of the former vibration is less than expected because in our setting the mono-and bilayer are simulated in a vacuum, without any interaction with typical substrates used in experimental work, such as silica or sapphire, that could affect the force constant of these modes.